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We have studied the space-time evolution of minimally viscous (^=0.08) QGP fluid, undergoing 
boost-invariant longitudinal motion and arbitrary transverse expansion. Relaxation equations for 
the shear stress tensor components, derived from the phenomenological Israel-Stewart's theory of 
dissipative relativistic fluid, are solved simultaneously with the energy-momentum conservation 
equations. Comparison of evolution of ideal and viscous fluid, both initialized under the similar 
conditions, e.g. same equilibration time, energy density and velocity profile, indicate that in viscous 
fluid, energy density or temperature of the fluid evolve slowly than in an ideal fluid. Transverse 
expansion is also more in viscous evolution. We have also studied particle production in viscous 
dynamics. Compared to ideal dynamics, in viscous dynamics, particle yield at high pr is increased. 
Elliptic flow on the other hand decreases. Minimally viscous QGP fluid, initialized at entropy density 
Sj n i=110 frn~ 3 at the initial time Ti=0.6 fm, if freeze-out at temperature Tf=130 MeV, explains the 
centrality dependence of pr spectra of identified particles. Experimental pr spectra of 7r~, K + and 
protons in 0-5%, 5-10%, 10-20%, 20-30%, 30-40% and 40-50% Au+Au collisions are well reproduced 
throughout the experimental pr range. This is in contrast to ideal dynamics, where, the spectra are 
reproduced only up to pr ~1.5 GeV. Minimally viscous QGP fluid, also explain the elliptic flow in 
mid-central (10-20%, 16-23%, 20-30%) collisions. The minimum bias elliptic flow is also explained. 
However, the model under-predict/over-predict the elliptic flow in very central/peripheral collisions. 

PACS numbers: 47.75.+f, 25.75.-q, 25.75.Ld 



I. INTRODUCTION 



Recent Au+Au collisions at the Relativistic Heavy Ion Collider (RHIC), at cm energy of GeV, give a 

strong indication that in central collisions, a hot dense matter is formed [H, 0, H, H[ . Whether the formed matter 
can be identified as the much sought after Quark-Gluon Plasma (QGP) as predicted in Lattice QCD simulations 
Q is presently debatable. However, experimentally observed elliptic flow in non-central Au+Au collisions gave a 
strong indication that a thermalized collective QCD matter is formed. Additionally, success of ideal fluid dynamics 
in explaining several experimental data e.g. transverse momentum spectra of identified particles, elliptic flow etc [||, 
together with the string theory motivated lower limit of shear viscosity 2 > 1 /47r 0, H| leading to a paradigm that in 
Au+Au collisions, a nearly perfect fluid is created. 

However, the paradigm of "perfect fluid" produced in Au+Au collisions at RHIC need to be clarified. It so happens 
that the ideal fluid dynamic models do have their shortcomings For example, experimentally, elliptic flow tends to 
saturate at large transverse momentum. The ideal fluid dynamics on the other hand predicts a continually increasing 
elliptic flow. The transverse momentum spectra of identified particles also starts to deviate from ideal fluid dynamics 
prediction beyond px ~ 1.5 GeV. Experimentally determined HBT radii are not reproduced in the ideal fluid dynamic 
models, the famous "HBT puzzle" [10j . Ideal fluid dynamics also works best in central collisions and gets poorer in 
more peripheral collisions. 

The shortcomings of ideal fluid dynamics possibly indicate greater importance of dissipative effects in the px ranges 
greater than 1.5 GeV or in more peripheral collisions. Indeed, ideal fluid is a concept, which is never realized in 
nature. As suggested in string theory motivated models, QGP viscosity could be small, rj/s > 1/4-7T, nevertheless it 
is non-zero. It is important to study the effect of viscosity, even if small, on space-time evolution of QGP fluid and 
quantify its effect. This requires a numerical implementation of relativistic dissipative fluid dynamics. Furthermore, if 
QGP fluid is formed in heavy ion collisions, it has to be characterized by measuring its transport coefficients, e.g. heat 
conductivity, bulk and shear viscosity. Theoretically, it is possible to obtain those transport coefficients in a kinetic 
theory model. However, in the present status of theory, the goal can not be achieved immediately, even more so for a 
strongly interacting QGP (sQGP). Alternatively, one can use the experimental data to obtain a "phenomenological" 
limit to the transport coefficients of sQGP. It will also require a numerical implementation of relativistic dissipative 
fluid dynamics. There is another incentive to study dissipative hydrodynamics. Ideal hydrodynamics depends on the 
assumption of local equilibrium. Dissipative hydrodynamics on the other hand do not depend on the assumption of 
local equilibrium. As will be discussed later, the fluid should not be far from an equilibrium state. The range of 
validity of dissipative hydrodynamics thus increases. Indeed, we can explore early times of fluid evolution better in a 
dissipative hydrodynamics. 

Theory of dissipative relativistic fluid has been formulated quite early. The original dissipative relativistic fluid 
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equations were given by Eckart [ll| and Landau and Lifshitz [12J. They are called 1st order theories. Formally, 
relativistic dissipative hydrodynamics are obtained from an expansion of entropy 4-current, in terms of dissipative 
fluxes. In 1st order theories, entropy 4-current contains terms linear in dissipative quantities. 1st order theory of 
dissipative hydrodynamics suffers from the problem of causality violation. Signal can travel faster than light. Causality 
violation is unwarranted in any theory, even more in a relativistic theory. The problem of causality violation is removed 
in the Israel-Stewart's 2nd order theory of dissipative fluid [131] ■ In 2nd order theory, expansion of entropy 4-current 
contains terms 2nd order in dissipative fluxes. However, these leads to complications that dissipative fluxes are no 
longer function of the state variables only. They become dynamic. The space of thermodynamic variables has to be 
extended to include the dissipative fluxes (e.g. heat conductivity, bulk and shear viscosity). 

Even though 2nd order theory is formulated some 30 years back, significant progress towards its numerical imple- 
mentation has only been made very recently 0, [H El O E El 0, El, H3, [H, HI [11 [M 113, 111, [11, [M 111 ill . 
At the Cyclotron Centre, Kolkata, we have developed a code " AZHYDRO-KOLKATA" to simulate the hydrodynamic 
evolution of QGP fluid including the effect of dissipation due to shear viscosity only. The code can simulate both 
the 1st order and the 2nd order dissipative hydrodynamics. Some results of AZHYDRO-KOLKATA for first order 
dissipative hydrodynamics have been published earlier [24], HI [H [2?| ■ In an earlier paper [27| , we have studied the 
2nd order dissipative hydrodynamics . The study was limited to the QGP phase only. More recently, the code was 
extended to include the phase transition [28]. In the present paper, we study, in detail, 2nd order dissipative hydro- 
dynamics evolution of QGP fluid including phase transition. We will also compare viscous hydrodynamics predictions 
with experimental data. It will be shown that the experimental pt spectra of identified particles and elliptic flow are 
much better explained in 2nd order dissipative hydrodynamics than in ideal hydrodynamics. 

The paper is organized as follows: In section [II] we shortly review the Israel-Stewart's [l3| phenomenological theory 
of relativistic dissipative fluid dynamics. We will limit our study to boost-invariant motion and consider dissipation 
only due to the shear viscosity. In section ULT1 we derive the relevant equations in 2+1 dimension. In section ITVl we 
discuss the equation of state, viscosity coefficient and initial conditions used in the present study. Hydrodynamic 
evolution gives the temporal evolution of thermodynamics variables, e.g. energy density or temperature, fluid velocity 
components, the shear stress tensor components. The information needs to be converted to particle spectra to compare 
against experiments. InfVl we present the relevant equations for computing invariant particle distribution in viscous 
dynamics. We test the numerical accuracy of the code "AZHYDRO-KOLKATA" in section fVTl Evolution of ideal and 
minimally viscous fluid is compared in section IVIII In section I Villi px spectra and elliptic flow in viscous dynamics 
are compared with the experimental data. Fluid evolution and subsequent particle production depend on the initial 
shear stress tensor and also on the relaxation time. Some preliminary study is described in section HXl The concluding 
section IXl summarizes our results. 



II. DISSIPATIVE FLUID DYNAMICS 



In this section, I briefly discuss the Israel-Stewarts phenomenological theory of dissipative hydrodynamics. More 
detailed exposition can be found in [l|. 

A simple fluid, in an arbitrary state, is fully specified by primary variables: particle current (iV M ), energy- momentum 
tensor (T^) and entropy current (S" M ) and a number of additional (unknown) variables. Primary variables satisfies 
the conservation laws; 

d»N» = 0, (2.1) 
d^ v = 0, (2.2) 

and the 2nd law of thermodynamics, 



dftS" > 0. (2.3) 

In relativistic fluid dynamics, one defines a time- like hydrodynamic 4- velocity, u M (normalized as u 2 = 1). One also 
define a projector, = g^" — u^u v , orthogonal to the 4-velocity (A M1/ ?i„ = 0). In equilibrium, an unique 4-velocity 
(w M ) exists such that the particle density (n), energy density (e) and the entropy density (s) can be obtained from, 

K= nuf± ( 2 - 4 ) 
T£ = eu"u v - P A^ (2.5) 

S2 q = su» (2.6) 
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An equilibrium state is assumed to be fully specified by 5-parameters, (n, e, m m ) or equivalently by the thermal 
potential, a = /i/T (/x being the chemical potential) and inverse 4-temperature, (3^ — u^/T. Given a equation of 
state, s — s(s,n), pressure p can be obtained from the generalized thermodynamic relation, 

S^pP-aN^+foTg (2.7) 

Using the Gibbs-Duhem relation, d(p(3^) — N£ q da — T^dpx, following relations can be established on the equilib- 
rium hyper-surface S eg (a,/3 M ), 

dS£ l = -adN*+0xdr% (2.8) 

In a non-equilibrium system, no 4- velocity can be found such that Eqs l2.4l2.5|2.6l remain valid. Tensor decomposition 
leads to additional terms, 

= Ng q + 5N>* = nu^ + V 1 (2.9) 
= [suV - pA**"] + IW + n"" 

+(iw + r/) (2.io) 

= Sg q + SS" = su^ 1 + ^ (2.11) 

The new terms describe a net flow of charge = A^ V N U , heat flow, = (e + p)/nV^ + q 11 (where g M is the 
heat flow vector), and entropy flow <K. II = —^A^ u T tJll/ — p is the bulk viscous pressure and t:^ = [\(A^ a A VT + 
A ucr A flT — ^A^" A aT ]T aT is the shear stress tensor. Hydrodynamic 4-velocity can be chosen to eliminate either V 1 
(the Eckart frame, u M is parallel to particle flow) or the heat flow (the Landau frame, u M is parallel to energy 
flow). In relativistic heavy ion collisions, central rapidity region is nearly baryon free and Landau's frame is more 
appropriate than the Eckart's frame. Dissipative flows are transverse to u^ 1 and additionally, shear stress tensor is 
traceless. Thus a non-equilibrium state require 1+3+5=9 additional quantities, the dissipative flows II, g M (or U M ) 
and 7r Miy . In kinetic theory, and T MI/ are the 1st and 2nd moment of the distribution function. Unless the function 
is known a-priori, two moments do not furnish enough information to enumerate the microscopic states required to 
determine S" M , and in an arbitrary non-equilibrium state, no relation exists between, N u ', T^ v and S^. Only in a 
state, close to an equilibrium one, such a relation can be established. Assuming that the equilibrium relation Eq |2.8l 
remains valid in a "near equilibrium state" also, the entropy current can be generalized as, 

S" = Sg q + dS" = p^ - aN» + (3 x T Xfl + (2.12) 

where Q» is an undetermined quantity in 2nd order in deviations, 6N>* = A" 1 - Ng q and ST**" = T^ v - Tg q . Detail 
form of Q M is constrained by the 2nd law d^S 11 > 0. With the help of conservation laws and Gibbs-Duhem relation, 
entropy production rate can be written as, 

8^ = -SN^d^a + ST^d^ + d^Q" (2.13) 

Choice of Q M leads to 1st order or 2nd order theories of dissipative hydrodynamics. In 1st order theories the 
simplest choice is made, = 0, entropy current contains terms up to 1st order in deviations, 5N 11 and ST 1 *". 
Entropy production rate can be written as, 

Td ll S» = TlX-q»X ll + ^ v X llv (2.14) 

where, X = -V.u; X^ = ¥f - u v d v u» and X» v = V<^u v> . 

The 2nd law, d^S^ > can be satisfied by postulating a linear relation between the dissipative flows and thermo- 
dynamic forces, 



n= -CO, (2.15) 
nT 2 

<f= -A— V'WT), (2.16) 
e +p 

= 2??V < V > (2.17) 
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where £, A and rj are the positive transport coefficients, bulk viscosity, heat conductivity and shear viscosity respec- 
tively. 

In 1st order theories, causality is violated. If, in a given fluid cell, at a certain time, thermodynamic forces vanish, 
corresponding dissipative fluxes also vanish instantly. Violation of causality is unwanted in any theory, even more 
so in relativistic theory. Causality violation of dissipative hydrodynamics is corrected in 2nd order theories [13(. In 
2nd order theories, entropy current contain terms up to 2nd order in the deviations, ^ 0. The most general 
containing terms up to 2nd order in deviations can be written as, 

& = -(AH 2 - 0iq v q v + /W^)^ - ^ + (2.18) 

As before, one can cast the entropy production rate (Td^S^) in the form of Eq |2.14l Neglecting the terms in- 
volving dissipative flows with gradients of equilibrium thermodynamic quantities (both are assumed to be small) and 
demanding that a linear relation exists between the dissipative flows and thermodynamic forces, following relaxation 
equations for the dissipative flows can be obtained, 



-C(0 + /3 DU) 
" nT 2 



-A 



P T 



(2.19) 
(2.20) 
(2.21) 



where D = u^d^is the convective time derivative. Unlike in the 1st order theories, in 2nd order theories, dynamical 
equations control the dissipative flows. Even if thermodynamic forces vanish, dissipative flows do not vanish instantly. 
It is important to mention that the parameters, a and (3\ are not connected to the actual state (iV M , T^ v ). The pressure 
P in Eq 12.121 is also not the "actual" thermodynamics pressure, i.e. not the work done in an isentropic expansion. 
Chemical potential a and 4-inverse temperature f3\ has meaning only for the equilibrium state. Their meaning need 
not be extended to non-equilibrium states also. However, it is possible to fit a fictitious "local equilibrium" state, 
point by point, such that pressure p in Eg 12. 121 can be identified with the thermodynamic pressure, at least up to 1st 
order. The conditions of fit fixes the underlying non-equilibrium phase-space distribution. 

It may be mentioned here that relaxation equations for the dissipative fluxes can also be derived in kinetic theory 
[3lH34j. Relaxation e quat ion from kinetic theory can contain additional terms, which are missed in the Israel-Stewart's 
theory. In [2(1 |2~l1. l29l l30j| . authors have used the following relaxation equation, derived from the kinetic theory, 



tt^ = 2?y [\7 < ^u" > - T n DTr^] - tvKtt"* + u^^Dux (2.22) 

which contain an additional term R = Du\ vis-a-vis Israel-Stewart's relaxation equation. The 

additional term ensures that throughout the evolution shear stress tensor remains traceless and transverse to fluid 
velocity. Israel-Stewart's theory is based on gradient expansion of entropy, gradients of equilibrium thermodynamic 
quantities (e, u) are assumed to be small. The term R does not contribute to entropy production and is missed 
in Eq l2.21l In the present work, we use Israel-Stewart's relaxation equation (neglecting the term R) to compute 
evolution of QGP fluid. As assumed in Israel-Stewart's theory, if gradients of thermodynamical quantities are small 
and since in principle, non-equilibrium effects are supposed to be small, both Du^ and it 111 ' small are small quantities 
and their product tt^Du^ can contribute only in 2nd order. It will be shown later that the contribution of the term 
R in hydrodynamic evolution of minimally viscous fluid is negligible. The evolution of energy density of the fluid is 
hardly changed whether or not the term R is included in the relaxation equation. However, for more viscous fluid, 
relaxation equation Eq |2.22l will be more appropriate than Israel-Stewart's equation 12.211 



III. (2+l)-DIMENSIONAL VISCOUS HYDRODYNAMICS WITH LONGITUDINAL BOOST 

INVARIANCE 



Complete dissipative hydrodynamics is a numerically challenging problem. It requires simultaneous solution of 14 
partial differential equations (5 conservation equations and 9 relaxation equations for dissipative flows). We reduce 
the problem to solution of 6 partial differential equations (3 conservation equations and 3 relaxation equations). In 
the following, we will study boost-invariant evolution of baryon free QGP fluid, including the dissipative effect due 
to shear viscosity only. Shear viscosity is the most important dissipative effect. For example, in a baryon free QGP, 
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heat conduction is zero and we can disregard Eq l2.20l Bulk viscosity is also zero for the QGP fluid (point particles) 
and Eq |2.f 91 can also be neglected. Shear pressure tensor has 5 independent components but the assumption of boost 
invariance reduces the number of independent components to three. For a baryon free fluid, we can also disregard 
the conservation equation Eq |2.U With the assumption of boost-invariance, energy-momentum conservation equation 
df 1 T m = become redundant. 

Heavy ion collisions are best described in (t, x, y, rf) coordinates, where r = \/t 2 — z 2 is the longitudinal proper 
time and r\ = h In |±| is the space-time rapidity. r± — (x, y) are the usual cartisan coordinate in the plane, transverse 
to the beam direction. With longitudinal boost-invariance the energy-momentum conservation equations T mn - n = 
yield 

d T f ^ + d X (T TT V x ) + dy(T TT Vy) = -( P + TV) (3.1) 
d T f ^ + d x (T TX V x ) + dy(T™V y ) = -d x (P + 0** - 7T TX V X ) - dy(^y - K^Vy) (3.2) 

d T f Ty + d x {f r y Vx ) + d y {f T y Vy ) = -d x {ir xy - n Ty v x ) - d y (p + n yy - n Ty v y ) (3.3) 

where A mn = tA"" 1 , p = rp, and v x = T TX /T TT , v y = T Ty /T TT . 

The components of the energy momentum tensor, including the shear pressure tensor are, 



T TT = (e+phl -p + tt tt (3.4) 
T*= {e+ P y f lv x +^ x (3.5) 
T Ty = {e + pYflvy + tt tv (3.6) 

We note that unlike in ideal fluid, in viscous fluid dynamics, conservation equations (see Eqs l3.Hl3.3p contain 
additional pressure gradients due to shear viscosity. Both T TX and T Ty components of energy-momentum tensor now 
evolve under additional pressure gradients. The rightmost term of Eq |3.1l also indicate that in viscous dynamics, 
longitudinal pressure is effectively reduced (note that the tt 1111 component is negative). Since pressure can not be 
negative, shear viscosity is limited by the condition, p + T 2 ir m > 0. The energy momentum conservation equations 
13. 1113.31 contain seven shear stress tensor components, tt xx , ir yv , ir xy , tt vv , ir TT , ir TX and 7r Ty . In 2nd order theory, the 
shear stress tensor is dynamical and evolves with time and the energy momentum conservation equations are to be 
solved together with the relaxation equations for the shear stress tensor, As mentioned earlier, with boost-invariance, 
only three stress tensor components are independent. We choose ir xx , tt vv and ir xy as the independent components. 
The dependent shear stress tensor components can be expressed in terms of the independent ones. The relaxation 
equations for the independent shear stress tensor components, in (r, x, y, rj) co-ordinates can be written as, 



d T n xx + v x d x n xx + v y d y Tr xx = — (n xx - 2 W XX ) (3.7) 

d T n yy + v x d x ir yy + v y d y ir yy = — (ji yy - 2 W yy ) (3.8) 

T?r7 

d T n xy + v x d x TT xy + v y d y TT xy = — (tt x « - 2ria xy ) (3.9) 

where r n is the relaxation time, r T = 2r]P2 (see Eq l2.21j) . The viscous pressure tensor relaxes on a time scale to 2r\ 

times the shear tensor = V'^tt"'. The xx, yy and xy components of the shear tensor can be written as 

a xx = -d x u x ~ u x Du x - -A xx 6 (3.10) 

3 

a yy = ~d y u y - u y Du y - -A yy 9 (3.11) 



o xy = — [d x u y - d y u x - u x Du y ~ u y Du x ] 



— A xy 9 (3.12) 
3 v ' 

Here Du^ = u^d^ is the convective time derivative and 9 = d^u^ is the expansion scalar. 

The dependent shear stress tensor components can easily be obtained from the independent ones. Using the 
properties that (i) tt^ is transverse to and (ii) tt^ is traceless, g^TTuu = 0, the dependent shear stress tensor 
components can be obtained as, 
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„-TX „xx I n . „xy 

7T = V X 7T + VyTT u 

n Ty = v x ir xy + v y ir vy 

tt tt = vItt xx + v 2 y TT vv + 2v x v y ir xy 

T 2 n r ' 7 ' = -(l-u2)7r xx -(l-^)7r ra 



(3.13) 
(3.14) 
(3.15) 



+2v x V v TT Xy 



(3.16) 



In causal dissipative hydrodynamics, energy momentum conservation equations [Eqs l3.1ll3.3| are solved simulta- 
neously with the relaxation equations [Eqs l3.7H3.9| . In ideal hydrodynamics, the energy-momentum conservation 
equations are closed with an equation of state and given initial conditions, e.g. energy density (e) and fluid velocity 
(v x and v y ) distributions at time n, the equation of motions can be integrated to obtain the variables at the next time 
step T i+ i. While for ideal hydrodynamics, this procedure works perfectly, viscous hydrodynamics poses a problem 
that shear stress tensor components contains time derivatives, d T "f±, d T u x , d T u x etc. Thus at the time step Tj one 
needs the still unknown time derivatives. Numerically, time derivatives at step Tj could be obtained if velocities at 
time step 7j and Ti+i are known. One possible way to circumvent the problem, is to use time derivatives of the 
previous step, i.e. use velocities at time step Tj-i and 7~i to calculate the derivatives at time step t% (23l |. It implicitly 
assume that the fluid velocity at the time step i is the average of the fluid velocity at time step i — 1 and i + 1. 



One of the most important inputs of a hydrodynamic model is the equation of state. Through this input, the 
macroscopic hydrodynamic models make contact with the microscopic world. In the present calculation we have used 
the equation of state, EOS-Q, developed in ref. Q. It is a two-phase equation of state. The hadronic phase of EOS-Q 
is modeled as a non-interacting gas of hadronic resonance. As the temperature is increased, larger and larger fraction 
of available energy goes into production of heavier and heavier resonances. This results into a soft equation of state, 
with small speed of sound, c 2 ps 0.15. With increasing temperature, the available volume is filled up with resonances 
and the hadronic states start to overlap, and microscopic degrees of freedom are changed from hadrons to deconfincd 
quarks and gluons. The QGP phase is modeled as that of a non-interacting quark (u,d and s) and gluons, confined 
by a bag pressure B. Corresponding equation of state, p = |e — |B is stiff with a speed of sound <? s = |. The two 
phases are matched by Maxwell construction at the critical temperature, T c = 164MeV, adjusting the Bag pressure 
£? 1 / 4 =230 MeV. As discussed in Q, ideal hydrodynamics with the equation of state EOS-Q, explain a large volume 
of RHIC Au+Au data. 



Shear viscosity coefficient (rj) of dense nuclear medium (QGP or resonance hadron gas) is quite uncertain. While in 
principle, shear viscosity coefficient can be calculated in kinetic theory, it is a complex problem, more so in the strong 
coupling limit. The problem is yet unsolved. Recently, using the ADS/CFT correspondence @, [1] it was conjectured 
that the ratio of shear viscosity over the entropy ratio is bounded form the lower side. In a strongly coupled gauze 
theory, N=4 SUSY YM, the shear viscosity can be evaluated as rj = ^N 2 T 3 and the entropy density is given by 

s = ?yN 2 T 3 . Thus in the strongly coupled field theory, 



IV. EQUATION OF STATE, VISCOSITY COEFFICIENT AND INITIAL CONDITIONS 



A. Equation of state 



B. Shear viscosity coefficient 




(4.1) 



ADS/CFT correspondence suggests a minimal viscosity for the strongly couple gauze theories. In the present paper, 
we have limited our study to the evolution of minimally viscous (rj/s = 1/An ps 0.08) QGP fluid. We assume that 
throughout the evolution, the ratio ?7/s=0.08 is maintained. 
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C. Relaxation time 



Relaxation time (t w ) is another important parameter in viscous hydrodynamics. As seen in Egs J3.7ti3.9l in the time 
scale tv shear stress tensor tt^ relaxes to 1st order values, ir^ = 2r\(j^ v . Instantaneous values of the shear stress 
tensor and the associated thermodynamical quantities, characterizing the fluid, will depend on the relaxation time. 
Again in principle, kinetic theory can be used to compute the relaxation time, but as with the viscosity, for a strongly 
coupled system, the complex problem is yet unsolved. In the Boltzmann gas limit, one can compute the relaxation 
time, 



4p Is 

where p is the pressure. In the present paper we have used this estimate for the relaxation time, However, we will 
also investigate the relaxation time dependence of the fluid evolution and subsequent particle production. 



D. Initial conditions 



Numerical integration of Eqs l3.Hl3~31 and 13. 7113. 91 require initial conditions, e.g. energy density e(x,y), velocity 
components v x (x, y) and v y (x,y), the shear stress tensor components, ir xx (x,y), ir vy (x,y) and tt xv (x, y) at the initial 
time Ti. A priori, the parameters arc unknown, only way to fix them is to confront the theory with experimental data. 

As discussed earlier, ideal hydrodynamics has been very successful in explaining a large volume of data in RHIC 
200AGeV Au+ Au collisions [6| . In [6| , the initial energy density profile was parameterized geometrically. At an impact 
parameter b, transverse distribution of wounded nucleons NwN(x,y,b) and of binary NN collisions NBc{x,y,b) 
are calculated in a Glauber model. A collision at impact parameter b is assumed to contain 25% hard scattering 
(proportional to number of binary collisions) and 75% soft scattering (proportional to number of wounded nucleons). 
Transverse energy density profile at impact parameter b is then obtained as, 



e{x, y, b) = £ (0.75 x N WN (x, y, b) + 0.25 x N BC (x, y, b)) (4.3) 

Assuming that at the initial time r^, fluid velocity is zero (v x (x,y) — v y (x,y) — 0), the parameter Eq and the 
initial equilibration time tj are fixed to reproduce the experimental transverse momentum distribution of pions in 
central Au+Au collisions. STAR and PHENIX data are fitted to obtain initial equilibrium time Ti — 0.6 fm and central 
entropy density of s = 110/m~ 3 . This corresponds to peak energy density of the fluid as 30 GeV/ fm 3 , or initial peak 
temperature of 350 MeV. 

In the present demonstrative calculation, we fix the initial time, energy density and fluid velocity, to the values 
obtained in ideal dynamics, i.e. at Ti=0.6 fm, Si n i=110fm~ 3 , v x = v y =0. In dissipative hydrodynamics, additionally, 
one need to specify the initial transverse profile for the three independent shear stress tensor components, ir xx , ir yy 
and TT xy . One choice could be to assume zero value initially, ir xx = ir yy — ir xy — 0. However, as it will be shown later, 
viscous effects are enhanced if initially shear stress tensor is non-zero. Since we are considering is boost-invariant 
evolution, one choice for non-zero shear stress tensor could be the boost-invariant value. We then assume that at the 
equilibration time Tj, the dissipative fluxes have attained their longitudinal boost-invariant values. The independent 
shear stress tensor components at the initial time is then obtained as, 



***(x,v)= (4.4) 

OTi 

n vv {x,y)= — , (4.5) 

OTi 

ir x v{x,y)= 0. (4.6) 

However, we do emphasize that the above choice is largely arbitrary. The correct values for initial 7r Miy can only be 
obtained by confronting the experimental data. 
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V. PARTICLE SPECTRA 

Solutions of Eqs l3~Tll3. 31 and !3. 7113.91 give the space-time evolution of the QGP fluid. The information about the time 
evolution of thermodynamic quantities e.g. energy density, fluid velocity and shear stress tensor, need to be converted 
into particle spectra to make contact with experiments. To do so, we use the well known Cooper- Frye prescription 35]. 
In Cooper-Frye prescription, the particle distribution is obtained by convoluting the one body distribution function 
over the freeze-out surface, 

where dH^ is the freeze-out hyper-surface and f(x,p) is the one-body distribution function. Now in ideal dynamics, 
the fluid is in local equilibrium and the one body distribution function is well approximated by the equilibrium 
distribution function, 

f( X ,p) = f°*(x )P ) = W( 1 -y—. (5.2) 

with inverse temperature f3 = 1/T and chemical potential /i. g is the degeneracy factor. In viscous dynamics on 
the other hand, the fluid is not in equilibrium and f(x,p) can not be approximated by the equilibrium distribution 
function f eq (x,p). In a highly non-equilibrium system, distribution function f(x,p) is unknown. If the system is 
slightly off-equilibrium, then it is possible to calculate correction to equilibrium distribution function due to (small) 
non-equilibrium effects. Slightly off-equilibrium distribution function can be approximated as, 

P e i(x,p)=ri(x 1 p)[l + cj)(x 1 p)], (5.3) 

4>{x,p) « 1 is the deviation from equilibrium distribution function f eq . With shear viscosity as the only dissipative 
forces, 4>(x,p) can be locally approximated by a quadratic function of 4-momentum, 

^{x,p)=E liV p" P v . (5.4) 
Without any loss of generality e^ u can be written as as, 

(5.5) 



2(e + P )T 2 ' 

completely specifying the non-equilibrium distribution function. As expected, correction factor increases with increas- 
ing viscosity. We also note that non-equilibrium correction depend quadratically on particle momentum. The effect 
of dissipation is more on large momentum particles. With the non-equilibrium distribution function thus specified, it 
is now possible to calculate the particle spectra from the freeze-out surface. In (r,x,y,r) s ) coordinates, the freeze-out 
surface is parameterized as, 

X/ 1 = (r/(x, y)coshrj s ,x,y,T f (x,y)smhr] s ), (5.6) 
and the normal vector on the hyper surface is, 

Qif dTf 

dS„ = (coshry s ,-- — ,-- — , - s'mh f] s )Tfdxdydri s (5.7) 
dx f dy f 

At the fluid position (r, x, y, ri s ) the particle 4-momenta are parameterized as, 

p' 1 = {m T cosh(rj s - Y),p x ,p y ,m T sinh(i] s - Y)) (5.8) 



The volume element p M dE M become, 
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p^dYi^ = (mTCOsh(r/ — Y) — pT-^TT^Tfdxdydr] (5-9) 
Equilibrium distribution function involve the term p ^ which can be evaluated as, 



p^u^ j(rriTCosh(r] — Y) — vt-Pt — m/t) 
T ~ T 

The non-equilibrium distribution function require the sum p^p^ir^, 



(5.10) 



VtxVv^ v = a\cosh 1 (r\ — Y) + a2COsh(rj — Y) + 03 (5-11) 

with 

01= m 2 T (Ti TT +t 2 tt w ) (5.12) 

a 2 = -2m T (p x TT TX +p y TT Ty ) (5.13) 

a 3 = p 2 x tt xx + p 2 y TT yy + 2p xPy ir xy - m 2 T T 2 ^ (5.14) 

Inserting all the relevant formulas in Eq l5.1l and integrating over spatial rapidity one obtains, 

dN _ dN e i + dN ne i ^ 
dyd 2 pT dyd 2 pT dyd 2 pr 

with, 



dN e * g 



dxdyTflrriTKiinP) - p T V T TfK (nf3)} (5.16) 



dyd 2 pT (2ir) 3 

dy^ ^WrJ dxd V T fi m T{^K,{nP) + ^K 2 (nf3) + (^ + a z )K x {np) + y#o («/?)} 



3 

-p T .V T T f {^-K 2 {n(3) + a 2 K x {n{3) + (y + a 3 )K {n{3)}] (5.17) 

where Kq, K\, K 2 and K3 are the modified Bessel functions. 

It must be mentioned that the non-equilibrium distribution function f neq (x,p) is obtained as a correction to the 
equilibrium distribution function f eq (x,p). It is implied that the non-equilibrium correction is small (f>(x,p) << 

dyd' 2 p T 

„. , ... TVi vnf in 

dyd 2 p T ' 



1. Consequently, non-equilibrium contribution to the invariant distribution must be small compared to the 

equilibrium contribution , dI ^ 2 - ■ The ratio, 



dyd A p T 

must be much less than 1. The condition limits the applicability of viscous dynamics. From Eqs l5.4l and l5.5[ we observe 
that the non-equilibrium correction depend quadratically on the transverse momentum and also on the thermodynamic 
variables (energy density or temperature and shear stress tensor) at the freeze-out surface. The px range over which 
viscous dynamics remains applicable thus depend on the freeze-out condition. 
We will also show results for elliptic flow v 2 . It is defined as, 

cos(2<£)# 

f 27r dN V ' 



Jo " dyd^pT ^ 



Expanding to the 1st order, elliptic flow as a function of transverse momentum can be obtained as, 
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v 2 (pt) = vl q (p T ) 



J pt dpT d(p J d(j)COs(2,(l)) p T dpTd(f) 



eq - j- j ^ r __________ 



\ J < pt dpT dip J ' pt dpT d<p / 

= v?( PT ) + v c 2 ° rr (pT), (5.20) 

where w^ 9 is the equilibrium elliptic flow and v^ rT is the correction due to the non-equilibrium effects. Non-equilibrium 
correction grow quadratically with px (see Ea l2.5p . on the other hand grow less than linearly with px- An 
important conclusion can be reached from Eq |5.20l asymptotically, viscous hydrodynamics do not predict saturation 
of elliptic flow. 



VI. STABILITY AND ACCURACY TEST FOR AZHYDRO-KOLKATA 



The energy-momentum conservation equations I3.H3.31 and the relaxation equations l3.8H3.10l are solved simultane- 
ously using the code, AZHYDRO-KOLKATA, developed at the Cyclotron Centre, Kolkata. AZHYDRO-KOLKATA 
was built upon the publicly available code " AZHYDRO" (3(| where energy- momentum conservation equations for the 
ideal fluid, undergoing arbitrary transverse expansion and boost-invariant longitudinal motion are solved using the 
SHASTA-FCT algorithm. AZHYDRO is extensively tested and results are published. We made extensive changes in 
the code to introduce dissipation due to shear viscosity. We have checked that AZHYDRO-KOLKATA reproduces 
the AZHYDRO results for dissipation free fluid. Equation of motion in viscous hydrodynamics dynamics can not be 
solved analytically. One can not compare numerical solutions with analytical results. Any numerical code has to live 
with certain uncertainty. One try to limit the uncertainty to an acceptable value. The general procedure for checking 
a numerical code is as follow: (i) results should be stable against change in integration step lengths, (ii) any symmetry 
in the system should not be destroyed due to numerical inaccuracy and (iii) numerical inaccuracy should not lead to 
unphysical maxima or minima. We test our code by checking whether these principles are violated or not. 

In FigUl we have shown the contour plot of energy density in x-y plane after an evolution of 8 fm in central Au+Au 
collision (b=0). The black, red and blue lines are for integration step lengths, dx = dy = 0.2/m, dr = 0.02/to, 
dx = dy = 0.2fm,dT = 0.01/m, and dx = dy — 0.1, dr = 0.01/m respectively. The evolution of energy density is 
stable against change in integration step lengths. Even by halving the integration step lengths, the constant energy 
density surfaces occur nearly at the same spatial position. For example, spatial position of the 0.1 GeV surface is 
uncertain by less than a few percent. In b=0 Au+Au collisions, the initial energy density distribution is symmetric 
with respect to x and y. Hydrodynamics evolution should not destroy the symmetry. Indeed, as seen in Fig[TJ the x-y 
symmetry is maintained in the evolution. We also note that the energy density falls smoothly with the radius. There 
is no unphysical excess or deficit of energy density in the x-y plane. This indicate that the numerical solution of the 
hydrodynamic equations do not lead to any unphysical maxima or minima. 

In Figl^l the constant temperature contours in x — r plane, at a fixed y=0 fm arc shown. Constant temperature 
contour plots also show that the AZHYDRO-KOLKATA results are very stable against change in integration step 
lengths. Early in the evolution, the temperature evolution does not indicate any dependence on the integration 
step lengths. Only at later times, that too at large distance from the center of the fluid, we find small change in 
temperature evolution. We are not showing but x and y components of the fluid velocity also show little dependence 
on the integration step lengths. 

Evolution of shear stress tensor components also does not show any significant dependence on the integration step 
lengths. In Figs|3]and|3J we have shown the contour plot of the tt xx and ir yv in x — r plane, in an impact parameter 
b=6.5 fm Au+Au collision. The ordinate is at a fixed value, y=0. Here again, evolution of tt xx and ir vv , with different 
integration step lengths, agrees with each other. It is expected also. Instantaneous values of energy density and 
velocity depend on the shear stress tensor. They are almost independent of integration step lengths. Naturally, shear 
stress tensor components do not show appreciable dependence on the integration step lengths. FigsHJH establish that 
the numerical solutions of the code, AZHYDRO-KOLKATA is stable against integration step lengths. It maintains 
the symmetry of the system and no unphysical maxima or minima arise due to numerical inaccuracy. 

We have performed another check. Early in the evolution, even when the fluid is undergoing arbitrary transverse 
and boost-invariant longitudinal expansion, at the center, the fluid is, least affected by the transverse motion. Early 
in the evolution, fluid at the center, will follow the equation of motion for boost-invariant longitudinal expansion. 
Causal viscous hydrodynamics for boost-invariant longitudinal expansion is well studied 16]. For equation of state 
p = l/3e = aT 4 and viscosity rj = bT 3 , a and b are constant, temperature evolution of the fluid can be obtained by 
solving two ordinary differential equations [l6l |. 
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dT T T- 3 <t> 



dr 3r 12ar 

d<f> 2aT$ $ / l 5 dT\ 8aT 4 



dr 3& 2 \t T dr 9r 



(6.1) 
(6.2) 



where $ = T 2 ir ,m . 

In FigO we have compared the AZHYDRO-KOLKATA solutions at the center (x=y=0) with the solutions of one 
dimensional scaling expansion. The dotted lines are temperature evolution of fluid undergoing one dimensional scaling 
expansion. They are obtained by solving Eqs l6.ll and l6.2l We have shown results for two initial temperature Ti=0.450 
GeV and 0.358 GeV, initial time is 7^=0.6 fm. The viscosity to entropy ratio is 7//s=0.08. The solid lines in FigfS] 
are solutions from AZHYDRO-KOLKATA. They depict the temperature evolution of the fluid at the center. For this 
particular plot, we have used Woods-Saxon type of distribution as the initial energy density distribution. Woods- 
Saxon type of distribution is rather flat . Centre of the fluid will be less affected by the transverse expansion in initial 
Woods-Saxon type of distribution. FigfS] indicate that at the center, evolution of temperature of the fluid, undergoing 
2+1 dimensional expansion, closely agree with the temperature evolution of the fluid in 0+1 dimension expansion. 
The difference is less than 3%. Lastly, in Fig[5]and[71 we have shown the temperature evolution of the fluid in x — t 
plane for ry/s=0. 08, 0.04, 0.02, 0.01 and 0, respectively. In Fig [S] and [7] the contours are drawn for fixed y=0 fm and 
y=5 fm respectively. Ideal hydrodynamics results are recovered as r\j s gradually reduces to zero. The results shown 
in Figs [Tj7] give us confidence the code "AZHYDRO-KOLKATA" correctly evaluates the viscous hydrodynamics. 

In Fig[TJ7|we have used Israel-Stewart's relaxation equation (£0 12. 21( 1 to compute the evolution of the fluid. As 
mentioned earlier, Eq |2.21l neglect the term R = — [u^ir" + u w tt^ x ]Du\ under the approximation that both n^" 
and Du^ are small and their product contribute only in 2nd order. It is important to check the accuracy of the 
approximation. We have simulated b=0 and b=6.5 fm Au+Au collision with and without the term R in the relaxation 
equation. Results are shown in Fig[5]and[5J where constant energy density contours in x — r plane are drawn at a fixed 
value of y=0. The black and blue lines are for fluid evolution with and without the term R in the relaxation equation. 
Both for b=0 and b=6.5 fm Au+Au collisions, at early time energy density evolve nearly identically irrespective of 
the term R, indicating that the term R contribute minimally to fluid evolution. Only at late time, Israel-Stewart's 
relaxation equation results in marginally slower evolution. The result is understood. Initially, even though ■n^ 1 ' is 
non-zero, the gradients of velocity is zero and the term R donot contribute. As the fluid evolve with time, velocity 
gradients grow but shear stress tensor -k^ v decreases. For minimally viscous fluid, velocity gradient do not grow to 
large value and the product ir^Du^ contribute minimally. 



VII. COMPARISON OF FLUID EVOLUTION IN IDEAL AND MINIMALLY VISCOUS FLUID 

In the following we compare the evolution of minimally viscous ( ^=0.08) fluid with ideal fluid evolution. For both 
the fluid, at the initial time Ti=0.6 fm, the peak entropy density is Si„i=110 fm~ 3 with the Glauber model transverse 
density profile (see Eq |4.3p and fluid velocity is zero, v x = v y =0. The shear stress tensor is assumed to attain the 
boost-invariant value, ir xx — ir yy — J^-, ir xy =0. In Fig. [101 evolution of energy density in a b=0 Au+Au collision, in 
viscous and in ideal dynamics is compared. In the left panels (a), (b) and (c), the constant energy density contours 
in viscous dynamics after evolution of 2.6 fm, 4.6 fm and 8.6 fm are shown. Their counterparts in ideal dynamics are 
shown in the right panels. Initially, energy density contours are identical both in ideal and in viscous simulations. As 
the fluid evolve, contour plots in ideal and viscous dynamics no longer remain identical. Evidently, energy density 
evolves slower in viscous dynamics. It is expected, viscosity opposes the expansion and cooling, the evolution is slowed 
down. 

To understand better the evolution of the viscous vis. a vis. ideal fluid, in Fig [Til and [T2l constant energy density 
contours in x — r plane are drawn. The black and blue lines are for ideal and viscous fluid respectively. In Fig llll 
contours are drawn at a fixed y=0 fm and in Fig ll2l at a fixed y=5 fm. Effect of viscosity is clearly seen, temperature 
evolution is slowed down. Viscous effect is comparatively large away from the center than near the center. It is 
understood. Velocity gradient generates viscosity. Velocity gradients are comparatively large away from the center. 
One also observes that transverse expansion seems to increase with viscosity. 

In viscous dynamics, fluid cools slower that in ideal dynamics. In Fig |13[ temperature evolution of viscous and 
ideal fluids at a fixed y=0 fm, but at different values of x=0,l,2,4, and 6 fm are shown. The black and blue lines 
corresponds to ideal and viscous fluid evolution. Viscosity enhances the life time of the QGP phase, mixed phase. For 
a fixed value of Tp, life-time of the hadronic phase is also increased. It is obvious that in an inhomogeneous medium, 
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effect of viscosity is not same everywhere. Viscous effects are comparativeiy less near the center of the fluid than at 
the periphery. 

Fluid velocity is assumed to be zero initially. With time fluid velocity grows. Evolution of fluid velocity in viscous 
and in ideal dynamics is compared in Figs ll4l and 1151 Contours of constant v x and v y in x-y plane, after 8 fm of 
evolution is shown in Fig [14] and Fig|T5j The black and blue lines are for ideal and viscous fluid respectively. Contour 
plots indicate that in the interior of the fluid, the fluid velocity grow comparatively faster in viscous dynamics than in 
ideal hydrodynamics. In Fig |161 in four panels, we have shown the temporal growth of the x-component of the fluid 
velocity (v x ). The black and blue lines are for ideal and viscous fluid. In the interior, fluid velocity starts to grow 
early. For example, at x = y = lfm, immediately after the evolution starts, the fluid velocity starts to grow. It grow 
for about 5 fm and then tends saturates. Away from the center, growth is more rapid. Also the fluid velocity grows 
to larger value than at the interior. Fig |16l indicate that in viscous dynamics, fluid velocity grow more rapidly than 
in ideal dynamics. While in viscous dynamics, the fluid cools slower, the fluid velocity grow faster. 

A. Evolution of shear stress tensor 

As the fluid evolve, the shear stress tensor 7r M1/ also evolve. Temporal evolution of shear stress tensor components 
ir xx and ir vv is shown in Fig[T71 Contour plot of ir xx and ir vv in x-y plane, after evolution of 2.6 fm, 4.6 fm and 8.6 
fm are drawn. We have initialized the shear stress tensor with boost-invariant values, ir xx = ir yy = 2t]/t. ir xx and 
ir vv rapidly decreases with time. Initially, both tt xx and ir vy have similar distribution, but after a few fm of evolution, 
the distribution of tt xx and tt vv starts to differ. One also observes that later in the evolution, tt xx or ir vv are stronger 
at the periphery than in the interior of the fluid. For example, after 4.6 fm of evolution, ir xx in the interior are less 
by a factor of 5 or more than the value at the periphery. Velocity gradients generate shear stress tensor and they 
are strongest on the periphery. It is also interesting to note the similarity between the stress tensor components 
tt xx and ir vv . As seen from the Eqs l3.10 l and 13. ill ir xx and ir vv are related by a transformation x — > y and y — > x. 
Fig] 1 71 correctly depict the relationship. We have not shown the evolution of the other independent component, ir xy . 
Initially, ir xy is zero. With time ir xy grow, but it do not grow to large value. ir xy remains much less than the other 
two independent components ir xx and tt vv . 

VIII. PT SPECTRA AND ELLIPTIC FLOW IN VISCOUS DYNAMICS 

In this section we will study the particle production in viscous dynamics. Space-time evolution of QGP fluid 
was solved using the computer code " AZHYDRO-KOLKATA" . Assuming that the hadrons freezes out at a fixed 
temperature Tp, we construct the freeze-out surface and use the Cooper- Frey prescription to calculate the px spectra 
and elliptic flow. We have ignored the resonance contribution to particle production. 

It was shown earlier that in AZHYDRO-KOLKATA, fluid evolution remains essentially unchanged if integration 
step lengths are halved. Nevertheless, fluid evolution does show small dependence on the integration step lengths. 
Let us first investigate the change in particle production when equation of motions are integrated with step lengths, 
(i) dx = dy=0.2 fm, dr=0.02 fm and (ii) dx = dy=0.1 fm, e?r=0.01 fm. With identical initial conditions, for Au+Au 
collisions at b=6.5 fm, we have solved minimally (j]/s=0.08) viscous hydrodynamics, with the two sets of integration 
step lengths, and computed pion yield from the freeze-out surface at Tp— 150 MeV. For the sake of comparison, we have 
also computed the yield in ideal dynamics, under similar conditions. In Fig |18f a) we have shown the px distribution of 
7r~. The solid and the dashed lines correspond to ideal and viscous dynamics respectively. Black and blue lines show 
the yield obtained with step lengths (i)dx = dy=0.2fm, dr=0.02fm and (ii)dx = ciy=0.1fm,ciT=0.01fm respectively. 
In ideal or viscous dynamics, transverse momentum distribution of tt~ do not show any significant dependence on 
integration step length. Whether the fluid evolution is calculated with step lengths (ljdx = dy = 0.2/to, dr — 0.02 fm 
or (ii)cir = dy = 0.1 fm, dr = 0.01/to, over nine order of magnitude, pion yield remain the same, The small variations 
the evolution with different integration step lengths, as observed in FigJTJll do not affect the pion transverse momentum 
distribution. 

In Fig |18f b). we have shown the elliptic flow in ideal and viscous dynamics. Elliptic flow is very sensitive observable. 
Unlike the px distribution of 7r~, both in ideal and viscous dynamics, elliptic flow of 7r~ do show small dependence 
on the integration step lengths. By halving the step lengths, elliptic flow increases. However, the increase is small, 
less than 10%. The results shown in Fig |18b again confirm that numerical evaluation of viscous dynamics, as done in 
AZHYDRO-KOLKATA is stable and reasonably accurate. 

FiglHl also illustrate the effect of viscosity on particle production. In viscous dynamics, px spectra is flattened. 
Compared to ideal dynamics, in viscous dynamics more particles are produced at large px- High px particles are 
produced from the high temperature phase and in viscous dynamics, the fluid remains in the higher temperature 
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phase for comparatively longer duration, enhancing the high px yield. The elliptic flow on the other hand decreases in 
viscous dynamics. Asymmetric pressure gradient generates elliptic flow. Apparently, in viscous dynamics, the shear 
stress tensor reduces the asymmetry in the pressure gradients. Indeed, appearance of additional pressure gradients 
in the energy-momentum conservation Eqs l3.2tf3.3l does suggest that asymmetry in pressure gradients is reduced in 
viscous dynamics. The increase of high px yield and decrease of elliptic flow in viscous dynamics appears to remedy 
the drawbacks of the ideal hydrodynamics. It is known that in ideal hydrodynamics, particle production at large px 
is under-predicted while elliptic flow is over-predicted. 

As mentioned earlier, viscous dynamics is applicable till the non-equilibrium contribution to equilibrium production 
is less than unity. In FigfT9l the ratio of non-equilibrium to equilibrium contribution to pion yield, for freeze-out 
temperatures TV=130,140,150 and 160 MeV, in a b=6.5 fm Au+Au collision, are shown. The ratio depends sensitively 
on the freeze-out temperature. For freeze-out temperature Tp=160,150 and 140 MeV, non-equilibrium contribution 
exceed the equilibrium contribution beyond px ~2.1, 2.5 and 3.5 GeV. For Tp=130 MeV, the applicability range 
exceed 5 GeV. The non-equilibrium correction depend on the thermodynamic variables on the freeze-out surface, 
e^v = 2 (e+p)T 2 7r ^ K /T 6 (assuming e — aT 4 ) , (see Eq l5.5|) . The correction depend on two opposing effects, as 
the fluid freezes out at lower and lower temperature, the correction increases due to 1/T 6 dependence. However, as Tp 
is lowered, the fluid evolve for longer time and 7r Miy decreases, the correction factor is decreased. 7r Miy decreases much 
faster than 1 /T 6 and the non-equilibrium correction is effectively reduced as the freeze-out temperature is lowered. 

It is also interesting to note from Fig ll9l that at low px, the non-equilibrium correction makes a negative contribution 
to the equilibrium production. As similar effect is seen in [l5| . where, in a blast-wave model, effect of shear viscosity 
on the freeze-out surface was studied. Effect of viscosity on low px and high px particles are quite different. At low 
px, viscous effect reduces the particle yield, while at high px, particle yield is increased. 



A. Comparison with experimental data 

1. Elliptic flow 

Let us now confront the minimally viscous hydrodynamics with the experimental data in Au+Au collisions at RHIC. 
Viscous hydrodynamics in 2+1 dimensions has a quite a few number of parameters, (i) the initial time (t~;), (ii) for 
a Glauber model type of transverse density profile, the peak initial energy density or equivalently the initial entropy 
density (si n i) (iii) the initial fluid velocity (v x and v y ) and (iv) the initial values for the three independent shear stress 
tensor components (ir xx , ir vy and ir xy ). Additionally, one require the freeze-out condition, the freeze-out temperature 
Tp. The shear viscosity and the relaxation time even though are calculable from kinetic theory, as stated earlier, 
the complex calculations are yet to be done and presently, they has to be treated as parameters again. Very large 
number of parameters effectively reduces the efficacy of viscous dynamics. Better fit obtained to the data in viscous 
dynamics may be attributed to the increased parameter space. In the following, we limit our study to minimally 
viscous dynamics, rj/s=0.08. We also fix the relaxation time to the Boltzmann gas estimate, tv = 6r//sT. Even then, 
the number of free parameters is large. We do not attempt to fit experimental data by varying all the parameters. As 
stated earlier, in ideal dynamics, if the QGP fluid is initialized at the initial time Ti—0.6 fm, to peak entropy density, 
Si n i = ll0 fm~ 3 , and if the initial fluid velocity is zero, for freeze-out temperature Tf=100 MeV, the RHIC data on 
identified particle px spectra and elliptic flow are reasonably well explained in a limited px range, px <1.5 GeV [||. 
To begin with we assume that the initial conditions of the minimally viscous QGP fluid remains same as in ideal 
hydrodynamics, e.g. (i) Ti=0.6 fm, (ii) s™=110 /m~ 3 (iii)v x = v y =0. We also assume that at the initial time Tj, 
the three independent shear stress tensor components have attained the boost-invariant values ir xx = ir yv = 2rj /3ri, 



7i-^=0. We then try to fit the PHENIX data [37| on elliptic flow in 16-23% Au+Au collisions, just by varying the 
freeze-out temperature (Tp). Being a ratio, elliptic flow is very sensitive to the details of the model. Hydrodynamics 
is better tested against the elliptic flow data than against the px spectra. 

Before we show the results of fit, it must be mentioned that, since entropy is generated in viscous dynamics, assuming 
same initial condition for the viscous fluid as in ideal fluid will lead to increase in multiplicity. Indeed, we have checked 
that in central collisions, compared to ideal dynamics, tt~ multiplicity is increased by ~ 28% in viscous dynamics. 
In other word, if ideal dynamics with r;=0.6 fm, Si„i=110 fm~ 3 can reproduce the experimental multiplicity data, 
viscous dynamics, initialized similarly, will over-predict the data. However, since, ideal hydrodynamics, generally 
under-predict the multiplicity in central collisions (e.g. see Fig. 5 in [38|h viscous dynamics will possibly give better 
description of multiplicity data in central collisions. 

The solid circles in FigfSD] are the PHENIX data on the transverse momentum dependence of the elliptic flow in 
16-23% Au+Au collisions. 16-23% Au+Au collisions roughly corresponds to b=6.5 fm collision. In Figf2TJl in four 
panels, for freeze-out temperatures, TV=130,140,150 and 160 MeV the elliptic flow in b=6.5 fm Au+Au collisions 
are shown. We have shown separately (i) the equilibrium contribution (the dash-dot line), (ii) the non-equilibrium 
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contribution (the dash-dot-dot line) and (iii) the total contribution (the solid line). The black arrows in the panels 
indicate the px range above which the non-equilibrium contribution exceed the equilibrium contribution and viscous 
dynamics become inapplicable. Several interesting points about elliptic flow in viscous hydrodynamics can be noted. 
Up to pt=2 GeV, the equilibrium elliptic flow (v^ 9 ) is hardly affected by changing Tp. Only beyond px—^ GeV, 
equilibrium elliptic flow increases with increasing Tp. Even then the increase is marginal, e.g. at px=3 GeV, 
increases by ~ 10% if Tp is increased from 130 to 160 MeV. One also note that the equilibrium flow agree with 
the PHENIX measurements up to px ~ 3 GeV. If non-equilibrium correction to elliptic flow is ignored, it would 
have appeared that the PHENIX data on elliptic flow in 16-23% centrality collisions is explained in the model for 
freeze-out temperature in the range 130-160 MeV. The non-equilibrium correction makes a negative contribution and 
reduces the elliptic flow. For Tp=140-160 MeV, the non-equilibrium contribution sufficiently reduces the elliptic flow 
and the total flow (equilibrium + correction) under-predict the PHENIX data. As the freeze-out temperature is 
lowered, magnitude of the non-equilibrium correction decreases. For Tp=130 MeV, non-equilibrium correction to the 
equilibrium flow is very small and PHENIX data are well explained. As seen in Fig[20l data up to pt «3.6 GeV are 
explained. For comparison, in Fig l20l we have shown the elliptic flow in ideal dynamics under identical conditions. In 
ideal dynamics, data are not explained beyond px=l-5 GeV. 

Even though, minimally viscous hydrodynamics describe the elliptic flow in 16-23% centrality Au+Au collisions, it 
appears that the experimental saturation of elliptic flow at large px is not explained. Beyond px=3.6 GeV, elliptic 
flow continue to increase in viscous dynamics, though the rate of increase is slowed down. Saturation of elliptic flow 
is possibly beyond viscous hydrodynamics. Indeed, above px—3 GeV, "jet" physics become important. Jets can 
influence the elliptic flow. Recently, in (39j . it was shown that elliptic flow is reduced in presence of jet quenching. 
Possibly viscous hydrodynamics, supplemented with jet physics will explain the saturation of elliptic flow at large px- 

While minimally viscous hydrodynamics, for freeze-out temperature Tp=130 MeV, give excellent description of 
elliptic flow in 16-23% centrality collisions, in other centrality ranges of collisions, the agreement is not so good. In 
FigHU PHENIX data 0, on the elliptic flow in 0-10%, 10-20%, 20-30% and 30-40% centrality Au+Au collisions, 
are shown. The solid lines in Figfjl] are viscous hydrodynamics predictions for elliptic flow in b=3.2, 5.7, 7.4 and 
8.7 fm Au+Au collisions. They roughly corresponds to 0-10%, 10-20%, 20-30% and 30-40% Au+Au collisions. The 
freeze-out temperature is Tp=130 MeV. Elliptic flow in 10-20% and 20-30% centrality collisions are well explained in 
the model. However the model under predict the elliptic flow in 0-10% centrality collisions and over predict the flow 
in 30-40% centrality collisions. At px=2-3 GeV, viscous dynamics predicts ~ 35% less elliptic flow in 0-10% centrality 
collisions and ~16% more elliptic flow in 30-40% centrality collisions. It appears that the minimally viscous dynamics 
explain the elliptic flow reasonably well only in mid-central collisions. In peripheral or in very central collisions, the 
predicted elliptic flow does not agree with the experiment, px dependence of the minimum bias V2 , on the other 
hand is excellently described in the model. In Fig |221 the filled circles show the px dependence of minimum bias V2, 
as measured by the STAR collaboration [4l|. Data extends upto px ~6 GeV. The solid line in FigfJUis the prediction 
from minimally viscous hydrodynamics. The experimental data are reproduced excellently. Even though the model 
is incapable of explaining the elliptic flow in very central or very peripheral collisions, it does explain the minimum 
bias v% . As shown here, minimally viscous hydrodynamics under-predict V2 in very central collisions and over-predict 
V2 in peripheral collisions. The two opposite effects are neutralized on averaging and minimum bias V2 is explained. 
However, we must remember that the present analysis is limited in the sense that the elliptic flow in 16-23% centrality 
collisions was fitted by varying only the freeze-out temperature. The other parameters of the model e.g. initial time, 
initial energy density, initial shear stress tensor were kept fixed. Viscous evolution as well as particle production 
depends on those parameters. It may be possible to explains the centrality dependence of elliptic flow, by tuning 
those parameters. Evidently, much more efforts will be needed to explain the centrality dependence of elliptic flow. 



2. Transverse momentum distribution 



Minimally viscous (r//s=0.08) QGP fluid, initialized to peak entropy density Sj n j=110 fm~ 3 at the initial time 
Ti=0.6 fm, if freezes-out at temperature Tp— 130 MeV, well reproduce the STAR data on the px dependence of 
minimum bias V2 as well as PHENIX data on 1)2 in 10-20%, 16-23%, 20-30% centrality Au+Au collisions. V2 in very 
central or very peripheral collisions, is not so well reproduced. Let us now confront the model predictions for px 
spectra of identified particles with the experimental data. In Figj23l PHENIX data 42j on the transverse momentum 
distribution of tt", in 0-5%,5-10%,10-20%,20-30%,30-40% and 40-50% centrality Au+Au collisions are shown. Solid 
lines in Fig |23l show the predictions from minimally viscous hydrodynamics, in Au+Au collisions at impact parameter 
b=2.3, 4.1, 5.7, 7.4, 8.7 and 9.9 fm. Roughly they corresponds to 0-5%,5-10%,10-20%,20-30%,30-40% and 40-50% 
centrality collisions. Model predictions are normalized by a factor N = 1.4. Minimally viscous hydrodynamics 
correctly reproduces the px spectra of 7r~ in all the centrality ranges of collisions. In contrast to ideal dynamics, 
where the transverse momentum spectra of 7r~ could be explained only up to px ~1.5 GeV, in minimally viscous 
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hydrodynamics, with the same initial conditions for the fluid, the spectra are explained right up to px—3 GeV. Let us 
remind that, apart from the normalizing factor of N=1.4, there is no free parameter in the model. Considering that 
resonance production is neglected, the normalizing factor N=1.4 seems reasonable. 

Minimally viscous hydrodynamics also correctly reproduce the px spectra of other identified particles e.g. K + and 
proton. In Fig l24l we have compared the minimally viscous hydrodynamic predictions for pt spectra of K + with the 
PHENIX data [42]. PHENIX data on K + production cover a lower px range, px ~ 2GeV . In all the centrality ranges 
of collisions, the PHENIX data on K + are excellently reproduced in viscous dynamics. Here again, a comparable fit 
could not be obtained in ideal dynamics. Viscous hydrodynamics also well reproduces the px spectra of protons. In 
Figj25] predictions for proton spectra are compared with the PHENIX experiment [13] • PHENIX collaboration could 
measure proton spectra over much extended px range, j»t=4.25 GeV. Minimally viscous hydrodynamics, explains the 
transverse momentum distribution of proton throughout the px range. Here again, comparable description to the 
data could not be obtained in ideal hydrodynamics. 

Minimally viscous hydrodynamics also reproduces the centrality dependence of mean px for 7r~, K + and protons, 
in central and mid-central Au+Au collisions. In Fig|26j the PHENIX data [42] on the centrality dependence of mean 
px for 7r~, K + and proton are shown. The solid lines are viscous dynamics predictions. For N par t >=100, < px >, for 
7r~, K + or protons are well reproduced. Only in peripheral collisions when N part < 100, viscous dynamics produces 
more mean px than in the experiment. Figures I20H26I clearly demonstrate the importance of viscosity. Even in 
minimally viscous dynamics, the viscous effect on particle production is considerable. Experimental data e.g. px 
spectra of identified particles and elliptic flow, are much better explained in minimally viscous dynamics than in ideal 
dynamics. 

IX. INITIAL SHEAR STRESS TENSOR AND RELAXATION TIME DEPENDENCE ON FLUID 

EVOLUTION AND PARTICLE PRODUCTION 

In the foregoing analysis, we have assumed that initially, the shear stress tensor (n^) has attained the boost- 
invariant value, tt xx — Ti yy — ^- 1 TT xy = 0. We have also used the Boltzmann gas estimate for the relaxation time 
(t,,-). Even though the choices seems to explain the PHENIX data on the centrality dependence of particle identified 
px spectra, the centrality dependence of the elliptic flow is not satisfactory. The fluid evolution and subsequent 
particle production depend on the initial 7r MI/ and also on tv. They can be tuned to obtain better fit to the data. 
Even though, we have not tuned them to fit the data, in this section, we study the dependence of fluid evolution and 
subsequent particle production on the initial value of ir' 11 ' and the relaxation time t w . 

A. Dependence of fluid evolution on initial shear stress tensor 

To investigate the dependence of the evolution, on the initial shear stress tensor tt^ , we have simulated b=6.5 fm 
Au+Au collisions with (i) boost-invariant values for the initial shear stress tensor, ir xx = ir vv — 2rj/3Ti, ir xy =0 and 
(ii) zero values for the initial shear stress tensor, ir xx — ir yv = tt xv ~Q. The other conditions of the fluid remaining the 
same, i.e. Ti=0.6 fm, Sj n j=110 /to~ 3 , v x — v y =0, and 7V = 6ry/sT. In Fig ]27l temporal evolution of the stress tensor 
component ir xx , for the two cases is compared. The blue and black lines are for initial ir xx = and ir xx — 2rj/3Ti 
respectively. The solid, dashed and medium dashed lines show the evolution at the fluid cell positions x=y=l,3 and 
5 fm respectively. When initially tt xx = 0, it grow rapidly to a maximum value around r=l fm. Beyond r=l fm, tt xx 
decreases and after 5-6 fm of evolution, reduces to very small value. On the other hand, when initially tt xx = 2r]/3ri, it 
decreases monotonically from the initial high value. A similar behavior is seen for the other independent component 
Tr vy . From Fig |27l one understands that, other conditions remaining the same, viscous effects on the evolution is 
enhanced if, initially, shear stress tensor has attained the boost-invariant value rather than zero value. It is also 
evident from Fig |281 In Fig |28i the temperature evolution in the two cases, is compared. Constant temperature 
contours, ini-r plane at a fixed y=0, are shown. Fluid evolution is slower with initial non-zero initial shear stress 
tensor. 



B. Initial shear stress tensor dependence of pr spectra and elliptic flow 

As shown in Fig |28l effects of viscosity on fluid evolution is enhanced if ~n^ v is initialized with the boost-invariant 
value rather than zero value. Since viscosity enhances high px production and reduces the elliptic flow, one expect 
more flattened px spectra and less elliptic flow in evolution with initially boost-invariant ir liv than in evolution with 
initially zero ■k ,1u . In Fig|29k and b, px spectra and elliptic flow in a b=6.5 fm Au+Au collisions, for (i) initial the 
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boost-invariant ir^ and (ii) initial zero tt^, are shown. Indeed, pr spectra is more flattened in evolution with initial 
boost-invariant 7r' Il/ (see Fig|2T)bL At Pt=3 GeV, a factor of 5 increase in 7r~ yield could be effected by changing the 
initial tt^ from zero to boost-invariant value. In Fig|29b. px dependence of elliptic flow is compared. Upto px ~ 1.5 
GeV, both the choice generate nearly same elliptic flow. Only beyond pt=L5 GeV, elliptic flow is less in evolution 
with initially boost-invariant it 111 ' than in evolution with initially zero 7r Ml/ . Compared to evolution with initially zero 
, at pr—i GeV, elliptic flow is reduced by ~ 20% in evolution with initially boost-invariant ir^ v . In Figl29b. we 
have also shown the PHENIX data on elliptic flow in 16-23% centrality Au+Au collisions. It appears that elliptic 
flow is better explained with initially boost-invariant . It is apparent that if initially shear stress tensor is zero, 
the data would require different freeze-out condition (or possibly different initial conditions). 

C. Dependence of fluid evolution on the relaxation time 

To investigate the effect of the relaxation time on fluid evolution, we have simulated Au+Au collision at b=6.5 fm, 
with three different values of the relaxation time, (i)r lr =-^, (ii) T~7r = fr and (m)T n =-]l. The other conditions of the 
fluid remaining the same. Fluid evolution, for the three relaxation time, is compared in Fig |30l In Fig |30l constant 
temperature contours in x — r plane at a fixed y=0 fm, are drawn. As the relaxation time increases, the fluid evolves 
slowly. Effect of relaxation time on the fluid evolution can be understood by examining the relaxation equation in 
one dimension, 

^=--(n-2 W ), (9.1) 

the solution of which is tt(t) ~ 1r\a + exp(— t/t^). Larger the relaxation time, more is the instantaneous values of 
the shear stress tensor. Naturally effect of viscosity is comparatively large with large relaxation time and the fluid 
evolves slowly with increasing relaxation time. 

D. Relaxation time dependence of pr spectra and elliptic flow 

In Fig l31b and b, we have compared the 7r _ pt spectra and elliptic flow in Au+Au collisions at b=6.5 fm, computed 
with three values of the relaxation time, 7V = 3r] / sT , 6r] / sT and 9r)/sT. The other condition of the fluid remaining 
the same. The pr spectra is flattened with increasing relaxation time (see Figl3Tk). Particle production at high pr 
is increased if the stress tensor takes longer to relax. The increase could be substantial. For example, at pr=3 GeV, 
changing relaxation time from ^ to |p, ir~ yield is increased nearly by a factor of 10. In Fig|3Tb. we have shown the 
elliptic flow i> 2 . Elliptic flow decreases as the relaxation time increases. For a factor of three increase in the relaxation 
time (^ = to |^), at pt ~ 3 GeV, elliptic flow decreases by ~ 15%. In Figl3"Tb. for comparison sake, we have 
shown the PHENIX data [37| on elliptic flow in 16-23% centrality collision. It appears that Boltzmann gas value 
7v = 6j]/sT, for the relaxation time best explain the PHENIX data on elliptic flow. 

X. SUMMARY AND CONCLUSIONS 

In Israel-Stewart's 2nd order theory of dissipative relativistic hydrodynamics, we have studied the evolution of QGP 
fluid in 2+1 dimensions. The "baryon free" QGP fluid (comprising u,d,s quarks and gluons) undergoes boost-invariant 
longitudinal motion and arbitrary transverse expansion. As the fluid expands its temperature decreases and when 
the temperature fall below a critical temperature T c =164 MeV, the fluid undergoes a 1st order phase transition to 
hadronic fluid (comprising of all the hadronic resonances). 

In 2nd order theory, in addition to the usual thermodynamic quantities e.g. energy density, pressure and hydro- 
dynamic velocity, dissipative flows are treated as extended thermodynamic variables. Relaxation equations for the 
dissipative flows are solved, simultaneously with the energy-momentum conservation equations. We have considered 
the most important dissipative effect, the shear viscosity. The bulk viscosity and heat conduction are neglected. 
ADS/CFT correspondence suggests a minimal viscosity for the QGP/hadronic fluid, ?7/s=0.08. We have limited our 
study to the evolution of minimally viscous QGP/hadronic fluid. 

In 2+1 dimensions, with boost-invariance, only three shear stress tensor components are independent. We choose 
it 11 , Tt yy and n xy as the independent components. The dependent components are obtained from the constraints on 
the shear stress tensor, e.g. tracelessness and transversality to fluid velocity. The formulation ensures that throughout 
the evolution, the constraints on the shear stress tensor are satisfied. The energy-momentum conservation equations 
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and the relaxation equations for the independent shear stress tensor components, are solved using a computer code 
" AZHYDRO-KOLKATA" , developed at the Cyclotron Centre, Kolkata. The code was tested extensively to ascertain 
its stability and accuracy. Explicit simulation of ideal and viscous fluids, initialized under similar conditions (e.g. 
same energy density, velocity profile at the equilibration time) confirms that the energy density or equivalently, the 
temperature of a viscous fluid, evolve slowly than its ideal counterpart. For a similar freeze-out condition, freeze-out 
surface is extended in viscous dynamics. Fluid velocity on the other hand, develops faster in viscous dynamics. 

As the fluid evolve, shear pressure tensor tt^ also evolve. We have assumed that at the initial time Tj, tt^ has 
attained the boost-invariant value. As the fluid evolve, the independent components n xx and n yy , which are non-zero 
initially, rapidly decrease and by 5-6 fm of evolution are reduced to very small values. The other independent shear 
stress tensor component ix xy is zero initially. It grows with time, but never grow to large value. Spatial distribution 
of tt xx and ir vv reveal an interesting feature of viscous dynamics. Even if, initially ir xx and ir vv have symmetric 
distribution, as the fluid evolve, pressure tensors quickly become asymmetric. However, in a central b=0, collision, 
we did not see any effect of the asymmetry in the energy density distribution. In a central collision, the two opposite 
asymmetry cancels each other. 

In viscous hydrodynamics fluid evolution depends considerably on the initial value of it^ 1 ' . If initially, instead of the 
boost-invariant value, tt^ is zero, effect of viscosity is much less on the evolution. For boost-invariant value, initially, 
tt xx — n yy — 2r)/ri, Tt xy =0. As the fluid evolve, ir xx and ir yy rapidly decrease. But in evolution with initial it 1 *" = 0, 
within a time scale of ~ 0.5 fm, n xx and Tr yy rapidly grow to reach a maximum and then decreases. The maximum 
reached is less than the corresponding boost-invariant value. Naturally, viscous effects are reduced. Viscous fluid 
evolution also depend on the relaxation time. Less the relaxation time, more is the effect of viscosity on evolution. 
For a given value of viscosity, viscous effect on the fluid evolution (and also on particle production) can be tuned by 
changing the relaxation time or the initial value of the shear stress tensor. 

We have also studied the effect of viscosity on particle production. Viscosity affects the particle production by 
(i) extending the freeze-out surface and (ii) by contributing a non-equilibrium correction d dyd » to the equilibrium 

particle d d d Ap T production. The non-equilibrium correction must be small, d dyd v « Jy%p T ' non_ec l um brium 
correction grow quadratically with pp and the condition d ^ d p q << dyd' 2 p T UIIU ts the pp range over which viscous 
hydrodynamics remains applicable. The non-equilibrium correction decreases as the fluid freeze-out at lower and 
lower temperature, pt range of applicability also increases accordingly. Both the effects (i) and (ii), enhances pt 
production and under similar conditions, compared to ideal dynamics, in viscous dynamics, particle yield at large pt 
is increased. The elliptic flow on the other hand decreases. Enhancement of pp production or reduction of elliptic 
flow in viscous dynamics can be understood in principle. In viscous dynamics, as the evolution is slowed down, the 
fluid spend more time in the higher temperature phase than in evolution in ideal dynamics, and the high pp yield 
is increased. In viscous dynamics, elliptic flow has two components, the equilibrium flow and the non-equilibrium 
correction to the equilibrium flow. The equilibrium flow itself is less than its ideal dynamics. The non-equilibrium 
correction makes a negative contribution to further reduce the flow. If we remember that asymmetry in the pressure 
gradients generate the elliptic flow, apparently, in viscous dynamics, asymmetry in pressure gradients are reduced. 
Indeed, appearance of additional pressure gradients in the energy-momentum conservation Eqs l3.2ll3~3l does suggest 
that asymmetry in the pressure gradients is reduced in viscous dynamics. 

Fitting experimental data in ideal hydrodynamics is a complex process, it is more complex in viscous hydrodynamics. 
As such ideal hydrodynamics has quite a few number of parameters, the initial time (r^), the (transverse) energy density 
and velocity profile at the initial time and a freeze-out condition (the freeze-out temperature Tp). In a remarkable 
effort, Kolb and Heinz [6J] fitted the RHIC data on particle identified pt spectra and elliptic flow in Au+Au collisions. 
For the transverse density profile they used the Glauber model calculation (see Eq |4.3p . It was seen that ideal 
QGP fluid, thermalised at Ti=0.6 fm, with peak entropy density Si„i=110 /to -3 and zero fluid velocity, v x =v y =0, if 
freeze-out at a temperature 1^=100 MeV, explains the data in a limited pp range, pp ~1.5 GeV. 

Parameter space is increased in viscous dynamics. As mentioned earlier, the three independent shear stress tensor 
components (ir xx , ir yy and ir xy ) and the viscosity coefficient (ry) and the relaxation time (r^) are to be treated as 
parameters. Since we consider only minimally viscous dynamics, the parameter r\ is fixed, 77/3=0. 08. We also fix 
the relaxation time to the Boltzmann gas value, tv = 6i]/sT. We also initialize the shear stress tensor 7r MI/ to 
boost-invariant value. To begin with we assume that in minimally viscous fluid, the initial condition of the fluid, 
i.e. equilibration time, initial entropy density, velocity, will not be changed much from the value obtained in ideal 
dynamics and initialize the fluid as in ideal hydrodynamics, i.e. Ti=0.6 fm, Si n i = 110 fm~ 3 , v x —v y =0. We only vary 
the freeze-out temperature (Tp) to fit PHENIX data [37| on the pp dependence of elliptic flow in 16-23% centrality 
Au+Au collisions. The PHENIX data on the pp dependence of elliptic flow in 16-23% Au+Au collisions can be well 
fitted if the hadronic fluid freeze-out at temperature 2f=130 MeV. Data up to pp ~3.6 GeV are then explained. 
With the same freeze-out temperature i.e. Tp=130 MeV, PHENIX data on the pp dependence of elliptic flow in 
mid central collisions e.g. 10-20%, 20-30% centrality collisions, are also well explained. The model also explains the 
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STAR data [4l| on the pr dependence of minimum bias V2- However, the elliptic flow, in very central or peripheral 
collisions, is not so well explained. While in very central collisions the model under predict the elliptic flow, in 
peripheral collisions, elliptic flow is over predicted. Apparently in very central collisions, PHENIX data on elliptic 
flow require less viscosity and in peripheral collisions data require more viscosity. Minimally viscous hydrodynamics 
with the parameters fixed from elliptic flow in 16-23% centrality Au+Au collisions, also explains the PHENIX data on 
pr spectra of identified particles, pr spectra of ir~ , K + and protons, in 0-10%, 10-20%, 20-30 — 5, 30-40% and 40-50% 
centrality collisions are explained more or less up to the experimental px range. It is a great achievement of viscous 
dynamics. As mentioned earlier, in ideal dynamics, data could be explained only up to pr ~1.5 GeV. The apparent 
contradiction i.e. minimally viscous hydrodynamics explains the centrality dependence of pt spectra of 7r~, K + and 
proton, but fails the explain the elliptic flow in very central or very peripheral collisions, indicate the sensitiveness of 
elliptic flow as compared to the p? spectra. Elliptic flow is a ratio and more sensitive to the details of the model than 
the pt spectra. Small difference between experiment and theoretical calculation, while unnoticed in the pt spectra 
gets amplified in the elliptic flow. Centrality dependence of the elliptic flow tests the model more stringently than the 
centrality dependence of the pr spectra. Minimally viscous hydrodynamics also reproduces the centrality dependence 
of mean transverse momentum of 7r~, K + and protons in central and mid-central collisions (up to A^ ar f=100). 

To conclude, present analysis demonstrated the importance of shear viscosity in RHIC Au+Au collisions. Even 
with the ADS/CFT lower bound on viscosity, the pr spectra and elliptic flow are modified to a great extent and the 
experimental data are better explained than in ideal hydrodynamics. However, it will not be proper to claim that 
RHIC has produced fluid as ideal as possible. A systematic analysis of RHIC data is required to check whether better, 
worse or equivalent description to the data are obtained with r/fs > 0.08. Only then one can claim about the QGP 
viscosity. 
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Figure Captions: 

Fig.l: (color online). Constant energy density contours in x-y plane after an evolution of 8fm. The black, red 
and blue lines are for integration step lengths: dx = dy=0.2 fm, e?r=0.02 fm, dx = dy=0.2 fm, dr=0.01 fm and 
dx = dy— 0.1 fm, <ir=0.01 fm. The viscosity coefficient is r]/s=0.08. 

Fig. 2: (color online), constant temperature contours in r — x plane. The black, red and blue lines are for integration 
step lengths, dx — dy=0.2 fm, dr=0.02 fm, dx = dy=0.2 fm, dr=0.01 fm and dx — dy— 0.1 fm, <iT=0.01 fm respectively. 

Fig. 3: (color online). Contour plot of n xx in x — r plane, for a fixed value of y=0 are shown. The viscosity coefficient 
is ?7/s=0.08. The black and blue lines are for integration step lengths, dx=dy=0.2 fm, dr=0.02 fm and dx=dy=0.1 
fm, dr=0.01 fm respectively. 

Fig. 4: (color online). Same as in Fig[3]but for the shear stress tensor tt vv . 

Fig. 5: The dotted lines depict the temperature evolution of viscous fluid in 0+1 dimension, for two initial tem- 
perature, Ti=.358 and .450 GeV. The solid lines are for the temperature evolution at the center of the fluid in 2+1 
dimensions. 

Fig. 6: Contour plot of temperature in x — r plane for fixed y=0. The five lines from top to bottom are for 
r]/s=0. 08, 0.04, 0.02, 0.01 and respectively. As the viscosity gradually reduces to zero, ideal hydrodynamics is recov- 
ered. 

Fig. 7: Same as in[6l but at y=5 fm. 

Fig. 8: (color online). Constant energy density contours in x — r plane at a fixed value of y=0 in a b=0 fm Au+Au 
collision. The black lines are obtained with relaxation equation Eg 12 . 221 containing the term R and the blue lines are 
obtained when R is neglected (see text). Effect of term R is not large. 

Fig. 9: (color online) same as in FigJ5]but for Au+Au collision at b=6.5 fm. 

Fig.l0:(color online). Evolution of energy density in minimally viscous (?y/s=0.08) and in ideal fluid evolution is 
compared. In panel a,b and c, constant energy density contours in x-y plane, at time r=2.6, 4.6 and 8.6 fm, in viscous 
evolution, are shown. Panels d,e and f show the energy density contours in ideal fluid evolution. 

Fig.ll: (color online). Constant energy density contours in x — t plane at a fixed y=0 fm. The black and blue lines 
corresponds to ideal and viscous fluid evolution. 

Fig. 12: (color online), same as in Fig[llJ but at y=5 fm. 

Fig. 13: (color online). Temporal evolution of fluid temperature at a fixed y=0 and different x-positions, x=0,l,2,4 
and 6 fm are shown. The black and blue lines are for ideal and viscous fluid evolution. Multiplicative factor used for 
different curves are shown in the figure. 

Fig.l4:(color online). Constant v x contours in x — y plane at t=8 fm. The black and blue lines correspond to ideal 
and minimally viscous fluid (7y/s=0.08) evolution. 

Fig.l5:(color online), same as in FigrjJ]but for v y . 

Fig: 16: (color online). Temporal evolution of the x-component of fluid velocity in ideal (the black lines) and in 
minimally viscous (the blue lines) dynamics is shown. 

Fig.l7:(color online). In left panels (a), (b) and (c) contour plot of the shear stress tensor component, tt xx in x — y 
plane at time r=2.6, 4.6 and 8.6 fm are shown. In right panels (d), (e) and (f) contours of constant ir yy are shown. 

Fig.l8:(color online). (a)Transverse momentum spectra for n~ . The blue and black lines are obtained when viscous 
dynamics is solved with integration step lengths (i)dx=dy=0.02,dr=0.02 and (ii) dx=dy=0.01,dr=0.01 respectively. 
The solid lines and dashed lines corresponds to ideal and minimally viscous (rj/s=0.08) fluid. ( b) Same as in (a) but 
for the elliptic flow. 

Fig. 19: The ratio of non-equilibrium contribution to equilibrium contribution to pion yield in b=6.5 fm Au+Au col- 
lisions. The solid, long dashed, medium dashed and short dashed lines are for freeze-out temperature 1^=130,140,150 
and 160 MeV. For applicability of viscous dynamics, the ratios must be much less than unity. 

Fig. 20: (color online). In four panels, px dependence of elliptic flow in b=6.5 fm Au+Au collisions, for freeze-out 
temperature Tf=160,150,140 and 130 MeV are shown. The dash-dot, dash-dot-dot and the solid lines are equilibrium 
elliptic flow, the non-equilibrium correction to the equilibrium flow and the total flow (equilibrium+ non-equilibrium 
correction), in minimally viscous hydrodynamics. The blues lines are elliptic flow in ideal hydrodynamics under similar 
conditions. The filled circles are the PHENIX data [37| on elliptic flow in 16-23% centrality Au+Au collisions. 

Fig.21:(color online). PHENIX data TO;] on the p T dependence of elliptic flow in 0-10%, 10-20%, 20-30% and 
30-40% Au+Au collisions are shown. The solid lines are predictions from minimally viscous hydrodynamics. 

Fig. 22: Filled circles are the STAR data [4l| on the pr dependence of minimum bias elliptic flow in Au+Au 
collisions. The solid line is the prediction from minimally viscous hydrodynamics. 

Fig.23:(color online). PHENIX data on p T spectra of tt" in 0-5%,5-10%,10-20%,20-30%,30-40% and 40-50% 
centrality Au+Au collisions are shown. The solid lines are predictions from minimally viscous hydrodynamics. 

Fig. 24: (color online). Same as in Fig|2"3l but for K + . 

Fig.25:(color online). Same as in Fig|2"31 but for proton. 
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Fig. 26: PHENIX data on the centrality dependence of average pr for n~ , K + and protons are shown. The solid 
lines are predictions from minimally viscous dynamics. 

Fig.27:(color online). Evolution of shear stress tensor tt xx (x, y = 0) for x=0,l,3 and 6 fm. The blue and black lines 
corresponds to fluid evolution with zero and non-zero (boost-invariant value) initial shear stress tensors. 

Fig.28:(color online). Dependence of fluid evolution on initial shear stress tensor. Constant temperature contours 
ill x — t plane at a fixed y=0 fm, in a Au+Au collision at b=6.5 fm are shown. The black lines correspond to initially 
zero shear stress tensor. The blue lines correspond to initially non-zero, boost-invariant shear stress tensors. 

Fig. 29: (a) Minimally viscous hydrodynamic predictions for the transverse momentum spectra of 7r~, in Au+Au 
collisions at impact parameter b=6.5 fm, for two values of initial shear stress tensors. The solid line corresponds to 
initial zero shear stress tensors Tr xx =Tr vv =TT xy —Q. The dashed line corresponds to initial boost invariant values for the 
shear stress tensors, / n xx =Tr yy =2r)/T i: tt xv =0. Pion yield decrease if initially shear stress tensors are zero, (b) same 
as in (a) but for the elli ptic flow for 7r~ . Elliptic flow increases if initially shear stress tensor is zero. The filled circles 
are the PHENIX data [37| on the transverse momentum dependence of elliptic flow in 13-26% centrality Au+Au 
collisions. 

Fig. 30: (color online). Dependence of fluid evolution on the relaxation time tv. The contours of constant temperature 
in x — t plane at a fixed y=0 fm, in Au+Au collision at b=6.5 fm are shown. The blue, red and black lines corresponds 
to relaxation time T w =3r)/sT, 6rj/sT and 9rj/sT. respectively. 

Fig.31:(a) Minimally viscous hydrodynamic predictions for the transverse momentum spectra of ir~ , in Au+Au 
collisions at impact parameter b=6.5 fm. The solid, dashed and short dashed lines corresponds to relaxation time, 
T n =3r]/sT, 6rj/sT and 9rj/sT. Pion yield increase as the relaxation time increase, (b) same as in (a) but for the 
ellip tic flow for ir~ . Elliptic flow decreases with increasing relaxation time. The filled circles are the PHENIX data 
[37j | on the transverse momentum dependence of elliptic flow in 13-26% centrality Au+Au collisions. 
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FIG. 1: (color online). Constant energy density contours in x-y plane after an evolution of 8fm. The black, red and blue lines 
are for integration step lengths: dx = dy— 0.2 fm, dr— 0.02 fm, dx = dy=0.2 fm, dr=0.01 fm and dx = dy=0.1 fm, dr=0.01 fm. 
The viscosity coefficient is rj/s=0.08. 
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FIG. 2: (color online), constant temperature contours in r — x plane. The black, red and blue lines are for integration step 
lengths, dx = dy=0.2 fm, dr—0.02 fm, dx — dy— 0.2 fm, dr=0.01 fm and dx = dy=0.1 fm, dr=0.01 fm respectively. 
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FIG. 3: (color online). Contour plot of ir xx in x — r plane, for a fixed value of y=0 are shown. The viscosity coefficient is 
?7/s=0.08. The black and blue lines are for integration step lengths, dx=dy=0.2 fm, dr=0.02 fm and dx=dy=0.1 fm, dr=0.01 
fm respectively. 
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FIG. 4: (color online). Same as in Fig|3]but for the shear stress tensor n 




FIG. 5: The dotted lines depict the temperature evolution of viscous fluid in 0+1 dimension, for two initial temperature, 
Ti=. 358 and .450 GeV. The solid lines are for the temperature evolution at the center of the fluid in 2+1 dimensions. 
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FIG. 6: Contour plot of temperature in x — r plane for fixed y=0. The five lines from top to bottom are for 
rj/s=0. 08, 0.04, 0.02, 0.01 and respectively. As the viscosity gradually reduces to zero, ideal hydrodynamics is recovered. 
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FIG. 7: Same as in [6] but at y=5 fm. 
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FIG. 8: (color online). Constant energy density contours in x — r plane at a fixed value of y=0 in a b=0 fm Au+Au collision. 
The black lines are obtained with relaxation equation Eq |2. 22! containing the term R and the blue lines are obtained when R 
is neglected (see text). Effect of term R is not large. 
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FIG. 9: (color online) same as in Fig|8]but for Au+Au collision at b=6.5 fm. 
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FIG. 10: (color online). Evolution ol energy density in minimally viscous (r?/s=0.08) and in ideal fluid evolution is compared. 
In panel a,b and c, constant energy density contours in x-y plane, at time t—2.6, 4.6 and 8.6 fm, in viscous evolution, are 
shown. Panels d,e and f show the energy density contours in ideal fluid evolution. 




FIG. 11: (color online). Constant energy density contours in x — t plane at a fixed y=0 fm. The black and blue lines corresponds 
to ideal and viscous fluid evolution. 




FIG. 12: (color online), same as in Fig llll but at y=5 



fm. 
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FIG. 13: (color online). Temporal evolution of fluid temperature at a fixed y=0 and different x-positions, x=0, 1,2,4 and 6 fm 
are shown. The black and blue lines are for ideal and viscous fluid evolution. Multiplicative factor used for different curves are 
shown in the figure. 




X (fm) 



FIG. 14: (color online). Constant v x contours in x — y plane at r~8 fm. The black and blue lines correspond to ideal and 
minimally viscous fluid (r;/s=0.08) evolution. 
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FIG. 16: (color online). Temporal evolution of the x-component of fluid velocity in ideal (the black lines) and in minimally 
viscous (the blue lines) dynamics is shown. 
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FIG. 17: (color online). In left panels (a), (b) and (c) contour plot of the shear stress tensor component, tt x 
time r=2.6, 4.6 and 8.6 fm are shown. In right panels (d), (e) and (f) contours of constant ir yy are shown. 
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FIG. 18: (color online). (a)Transverse momentum spectra for ir~ . The blue and black lines are obtained when viscous dynamics 
is solved with integration step lengths (i)dx=dy=0.02,dr=0.02 and (ii) dx=dy=0.01,dr=0.01 respectively. The solid lines and 
dashed lines corresponds to ideal and minimally viscous (r;/s=0.08) fluid. ( b) Same as in (a) but for the elliptic flow. 




FIG. 19: The ratio of non-equilibrium contribution to equilibrium contribution to pion yield in b=6.5 fm Au+Au collisions. 
The solid, long dashed, medium dashed and short dashed lines are for freeze-out temperature Tf = 130,140,150 and 160 MeV. 
For applicability of viscous dynamics, the ratios must be much less than unity. 
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FIG. 20: (color online). In four panels, pt dependence of elliptic flow in b=6.5 fm Au+Au collisions, for freeze-out temperature 
Tf=160,150,140 and 130 MeV are shown. The dash-dot, dash-dot-dot and the solid lines are equilibrium elliptic flow, the 
non-equilibrium correction to the equilibrium flow and the total flow (equilibrium+ non-equilibrium correction), in minimally 
viscous hydrodynamics. The blues lines are elliptic flow in ideal hydrodynamics under similar conditions. The filled circles are 
the PHENIX data [37[ on elliptic flow in 16-23% centrality Au+Au collisions. 
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FIG. 21: (color online). PHENIX data [lj on thep T dependence of elliptic flow in 0-10%, 10-20%, 20-30% and 30-40% Au+Au 
collisions are shown. The solid lines are predictions from minimally viscous hydrodynamics. 
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FIG. 22: Filled circles are the STAR data [411 ] on the pr dependence of minimum bias elliptic flow in Au+Au collisions. The 
solid line is the prediction from minimally viscous hydrodynamics. 
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FIG. 23: (color online). PHENIX data [42j on p T spectra of tt" in 0-5%,5-10%,10-20%,20-30%,30-40% and 40-50% centrality 
Au+Au collisions are shown. The solid lines are predictions from minimally viscous hydrodynamics. 



35 



10 3 i- 
10 2 

io- 

10° ! 

10" 2 i 

10 3 1 

10- i 

10 b i 

10 6 

10'' i 

10" 8 : - 




K+ in Au+Au@200 GeV(PHENIX) 
S ini =110fm- 3 ,Ti=0.6fm,T F =130MeV,ri/s=0.08 

5% 

10%(x.1) 




10 

20-30%(x.001) 
30-40%(x.0001) 
40-50%(x.00001) 

~> 1 1 i <~ 

1 



PT (GeV) 

FIG. 24: (color online). Same as in Fig[23]but for K + . 
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FIG. 25: (color online). Same as in Fig |23l but for proton. 
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FIG. 26: PHENIX data on the centrality dependence of average pr for tv~ , K + and protons are shown. The solid lines are 
predictions from minimally viscous dynamics. 
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FIG. 27: (color online). Evolution of shear stress tensor n xx (x,y — 0) for x=0,l,3 and 6 fm. The blue and black lines 
corresponds to fluid evolution with zero and non-zero (boost-invariant value) initial shear stress tensors. 
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FIG. 28: (color online). Dependence of fluid evolution on initial shear stress tensor. Constant temperature contours in x — r 
plane at a fixed y=0 fm, in a Au+Au collision at b=6.5 fm are shown. The black lines correspond to initially zero shear stress 
tensor. The blue lines correspond to initially non-zero, boost-invariant shear stress tensors. 
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FIG. 29: (a) Minimally viscous hydrodynamic predictions for the transverse momentum spectra of n~ , in Au+Au collisions 
at impact parameter b=6.5 fm, for two values of initial shear stress tensors. The solid line corresponds to initial zero shear 
stress tensors n xx =TT yy —n xy —0. The dashed line corresponds to initial boost invariant values for the shear stress tensors, 
tt xx =tt vv —2ri / Ti , tt xv —0. Pion yield decrease if initially shear stress tensors are zero, (b) same as in (a) but for the elliptic 
flow for 7r~. Elliptic flow increases if initially shear stress tensor is zero. The filled circles are the PHENIX data [37j on the 
transverse momentum dependence of elliptic flow in 13-26% centrality Au+Au collisions. 
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FIG. 30: (color online). Dependence of fluid evolution on the relaxation time 7V. The contours of constant temperature in x — r 
plane at a fixed y=0 fm, in Au+Au collision at b=6.5 fm are shown. The blue, red and black lines corresponds to relaxation 
time TTr=3ri/sT, &r)/sT and 9rj/sT. respectively. 
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FIG. 31: (a) Minimally viscous hydrodynamic predictions for the transverse momentum spectra of ir~ , in Au+Au collisions at 
impact parameter b=6.5 fm. The solid, dashed and short dashed lines corresponds to relaxation time, r^^Sr] / sT , G-q/sT and 
9rj/sT. Pion yield increase as the relaxation time increase, (b) same as in (a) but for the elliptic flow for 7r~. Elliptic flow 
decreases with increasing relaxation time. The filled circles are the PHENIX data [37| on the transverse momentum dependence 
of elliptic flow in 13-26% centrality Au+Au collisions. 



